load("partyleaders.RData")

names(leader)

###CODE TO RESET BASELINE CATEGORIES

reset.base<-function(x,base){
  if(!is.factor(x)){x<-factor(x)}
  x.lev<-levels(x)
  if(sum(1*(x.lev==base))==0){return(x)}
  ind<-which(x.lev==base)
  if(ind==1){return(x)}
  L<-length(x.lev)
  if(ind==L){return(factor(x,levels=x.lev[c(L,1:(L-1))]))}
  return(factor(x,levels=x.lev[c(ind,1:(ind-1),(ind+1):L)]))
}

#Make social democrats the baseline party family

leader$rbgroup<-reset.base(leader$groups, "30")
leader$groups<-leader$rbgroup


###MODEL 1-Time to First Female Leader (FFL)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

library(xtable)
xtable(summary(model1),digits=3)
round(exp(summary(model1)$coefficients[, 1]),3)

###
#MODEL 2-Party Leader Death 
####

summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)
round(exp(summary(deathmodellinear)$coefficients[, 1]),3)


###
##FUNCTIONS FOR INTERPRETATION WITH PREDICTED VALUES 
###

#loading mvtnorm for plotting
library(mvtnorm)

#computes empirical hpd
emp.hpd<-function(x, conf){
    conf <- min(conf, 1 - conf)
    n <- length(x)
    nn <- round(n * conf)
    x <- sort(x)
    xx <- x[(n - nn + 1):n] - x[1:nn]
    m <- min(xx)
    nnn <- which(xx == m)[1]
    return(c(x[nnn], x[n - nn + nnn]))
}

pred.boot.hazard<-function(model,x.mat,x.values,conf.level,unlink,num.draws,add=FALSE,to.unlink=TRUE,ylim,col,...){
    #model is an lm object
    #x.mat is a matrix of covariate values (gotten from setx for example)
    #where one value changes and all others are held constant
    #x.values is the vector of values for the covariate that is changing
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.95}
    if(missing(unlink)){unlink<-function(x){exp(x)}}
    if(to.unlink==FALSE){unlink<-function(x){x}}
    if(missing(num.draws)){num.draws<-10000}
   	beta.hat<-coef(model)
	beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred<-unlink(beta%*%t(x.mat))
    hpd.pred<-apply(pred,2,emp.hpd,conf=conf.level)
    med<-apply(pred,2,median)
    ymin<-min(cbind(med,t(hpd.pred)))
    ymax<-max(cbind(med,t(hpd.pred)))
    if(missing(ylim)){ylim=c(ymin,1)}
    if(missing(col)){col="black"}
    if(add==FALSE){
	    plot(x.values,med,ylim=ylim,col=col,...)
	}else{
		points(x.values,med,ylim=ylim,col=col,...)
	}
    points(x.values,hpd.pred[1,],type="l",lty=2,col=col)
    points(x.values,hpd.pred[2,],type="l",lty=2,col=col)
    return(list(x.values=x.values,median.pred=med,hpd.pred=hpd.pred))
}

diff.boot.hazard<-function(model,x.1,x.2,conf.level,unlink,num.draws){
    #model is an lm object
    #x.1 is one vector of covariate values (gotten from setx for example)
    #x.2 is one vector of covariate values (gotten from setx for example)
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.95}
    if(missing(unlink)){unlink<-function(x){1-exp(-exp(x))}}
    if(missing(num.draws)){num.draws<-100000}
   	beta.hat<-coef(model)
   	ind.na<-is.na(beta.hat)
   	beta.hat<-beta.hat[!ind.na]
   	x.1<-x.1[!ind.na]
   	x.2<-x.2[!ind.na]
    beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred.1<-unlink(beta%*%x.1)
    pred.2<-unlink(beta%*%x.2)
    median.1<-median(pred.1)
    median.2<-median(pred.2)
    d<-pred.2-pred.1 #difference of predicted values
    hpd.diff<-emp.hpd(d,conf.level)
    med.diff<-median(d)
    return(list(median.1=median.1,median.2=median.2,median.diff=med.diff,hpd.diff=hpd.diff))
}

db.2<-function(x.1,model,x.2,conf.level,unlink,num.draws){
    #model is an lm object
    #x.1 is one vector of covariate values (gotten from setx for example)
    #x.2 is one vector of covariate values (gotten from setx for example)
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.95}
    if(missing(unlink)){unlink<-function(x){exp(x)}}
    if(missing(num.draws)){num.draws<-100000}
   	beta.hat<-coef(model)
   	ind.na<-is.na(beta.hat)
   	beta.hat<-beta.hat[!ind.na]
   	x.1<-x.1[!ind.na]
   	x.2<-x.2[!ind.na]
    beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred.1<-unlink(beta%*%x.1)
    pred.2<-unlink(beta%*%x.2)
    median.1<-median(pred.1)
    median.2<-median(pred.2)
    d<-pred.1/pred.2 #cox hazard ratio
    hpd.diff<-emp.hpd(d,conf.level)
    med.diff<-median(d)
quantile.diff<-quantile(d,c((1-conf.level)/2,1-(1-conf.level)/2))
    return(list(median.1=median.1,median.2=median.2,median.diff=med.diff,hpd.diff=hpd.diff,quantile.diff=quantile.diff))
}

build.x.mat<-function(model,data.x.mat){
tt<-terms(model)
Terms<-delete.response(tt)
 m <- model.frame(Terms, data=data.x.mat, na.action = na.omit, 
            xlev = model$xlevels)
X <- model.matrix(Terms, m, contrasts.arg = model$contrasts)
return(X)
}


#####
##GENERATING PLOTS USING PREDICTED VALUES FOR MODEL 1-FIRST FEMALE LEADER
######

###Setting Values
quantile(leader[rownames(model.matrix(model1)),"seatdiff"],.1)
quantile(leader[rownames(model.matrix(model1)),"seatdiff"],.9)
summary(leader[rownames(model.matrix(model1)),"perfemmp2"])
x.values<-seq(-.5,.3,by=.1)#seatdiff values
Time<-20
opp=0
nogov=0
gov=1
seatdiff=x.values
seatdiff.mean=-0.0035
seatdiff.losing=-0.1
seatdiff.winning=0.1
seatdiffnochange=0
groups=30
groups4=30
perfemmp2.mean=16.590
perfemmp2.first= 3.846
perfemmp2.third= 27.300
seatshare=0.2193
minorparty=0
majorparty=1
Country="Finland"
#Country="Denmark"
#Country<-factor(Country,levels=levels(priya$Country))
FoundedAfter1980=0


##################
###Gov & Maj Party
#########

#Minor Party out of Gov
data.x.mat.oppmin<-data.frame(Time=1:48,seatdiff=seatdiff.mean,groups=groups4,perfemmp2=perfemmp2.mean,seatshare=seatshare,Country=Country,majorparty=minorparty,gov=nogov, FoundedAfter1980=FoundedAfter1980)

x.mat.model1<-build.x.mat(model=model1,data.x.mat.oppmin) #lowest
yoppmin<-cumprod(1-predict(model1,data.x.mat.oppmin,type="response"))
plot(cumprod(1-predict(model1,data.x.mat.oppmin,type="response")),type="l",lty=2,lwd=2,xlab="Time",ylab="Survival Probability",col="gray8") #red

#Minor Party in Gov
data.x.mat.govmin<-data.frame(Time=1:48,seatdiff=seatdiff.mean,groups=groups4,perfemmp2=perfemmp2.mean,seatshare=seatshare,Country=Country,majorparty=minorparty,gov=gov, FoundedAfter1980=FoundedAfter1980)


x.mat.model1<-build.x.mat(model=model1,data.x.mat.govmin) 

ygovmin<-cumprod(1-predict(model1,data.x.mat.govmin,type="response"))
points(cumprod(1-predict(model1,data.x.mat.govmin,type="response")),type="l",lty=1,lwd=2, col=1) #black

#Major Party in Gov
data.x.mat.govmaj<-data.frame(Time=1:48,seatdiff=seatdiff.mean,groups=groups4,perfemmp2=perfemmp2.mean,seatshare=seatshare,Country=Country,majorparty=majorparty,gov=gov, FoundedAfter1980=FoundedAfter1980)


x.mat.model1<-build.x.mat(model=model1, data.x.mat.govmaj)

ygovmaj<-cumprod(1-predict(model1,data.x.mat.govmaj,type="response"))
points(cumprod(1-predict(model1,data.x.mat.govmaj,type="response")),type="l",lty=3,lwd=2, col="gray60") #green

#Major Party out of Gov
data.x.mat.oppmaj<-data.frame(Time=1:48,seatdiff=seatdiff.mean,groups=groups4,perfemmp2=perfemmp2.mean,seatshare=seatshare,Country=Country,majorparty=majorparty,gov=nogov,FoundedAfter1980=FoundedAfter1980)


x.mat.model1<-build.x.mat(model=model1,data.x.mat.oppmin) #same as

yoppmaj<-cumprod(1-predict(model1,data.x.mat.oppmaj,type="response"))
points(cumprod(1-predict(model1,data.x.mat.oppmaj,type="response")),type="l",lty=4, lwd=2, col="gray48")#blue

legend(1, .3, c("Major Party, In Gov.", "Major Party, Out of Gov.", "Minor Party, In Gov.", "Minor Party, Out of Gov."), 
                col=c("gray60","gray48", "black","gray8"), lty=c(3,4,1,2 ), lwd=c(2,2,2,2))
                
########
#Survival Probs for Different Years (Gov/Opp and Major/Minor)
 yoppmin[(c(1,10,20,30))]
 ygovmaj[(c(1,10,20,30))]
 (ygovmaj-yoppmin)[(c(1,10,20,30))]
                
 
                
##################
###Seat Share
#########

###
##Basic Interpretation
####

exp(-0.092*10)# Interpretation for 10% change in seat share; exp of coeffiicent estimate times percentage gain in seats; used in inital analysis of first model

sd.vec=seq(-.20,.20,by=.01)
x.mat.bb<-data.frame(Time=1,seatdiff=sd.vec,groups=groups4,perfemmp2=perfemmp2.mean,seatshare=seatshare,Country=Country,majorparty= majorparty,gov=nogov,FoundedAfter1980=FoundedAfter1980)
x.mat.bb<-build.x.mat(model1,x.mat.bb)

x.bb<-apply(x.mat.bb,1,db.2,model=model1,x.2=x.mat.bb[sd.vec==0,],unlink=function(x){exp(x)},conf.level=.9, num.draws=100000)
x.bb<-matrix(unlist(x.bb),nrow=7)
plot(x.bb[3,]~I(sd.vec*100),type="l",col="black",lty=1,lwd=2,ylim=c(min(x.bb[4:5,]),max(x.bb[4:5,])),xlab="Percentage Difference in Seat Share",ylab="Cox Hazard Ratio")
#In tex or word, you would want to call this plot 'plot of the cox hazard ratio (y-axis) between the change in seat share (x-axis) and change in seat share equal to XXXX - the left most value of percentage difference in seat share'
points(x.bb[6,]~I(sd.vec*100),type="l",col="black",lty=2)
points(x.bb[7,]~I(sd.vec*100),type="l",col="black",lty=2)


#####
##GENERATING  COX HAZARD RATIOS USING PREDICTED VALUES FOR MODEL 1-FIRST FEMALE LEADER
###Use db.2 to get cox ratios
##Use hpds for CIs

#example using diff.boot
Time<-c(20,20)
nogov<-c(0,0)
gov<-c(1,1)
#majorparty<-c(0,1)
minorparty<-c(0,0)
majorparty<-c(1,1)
#majorparty<-c(1,0)
#gov<-c(1,0)
#majorparty<-c(1,1)
#seatdiff<-c(.38,-.57) #max and min
#seatdiff<-c(0.0261,-0.0336) #quartile
#seatdiff<-c(0.07395498,-0.07565408) #decile
seatdiff.10<-c(0.10,-.10) 
seatdifften<-c(0.10,0)
#seatdiff<-c(0.15,-.15) 
seatdiff<-c(0.00,0.00)
groups4<-c(30,30)
groups.switch<-c(30,40)
#groups4<-factor(groups,levels=levels(priya$groups))
perfemmp2<-c(16.590,16.590)
seatshare<-c(0.2193,0.2193)
#seatshare<-c(0.18340,0.18340)
age<-c(50,50)
#age<-c(86,86)
Country="Finland"
#Country<-factor(Country,levels=levels(priya$Country))
FoundedAfter1980<-c(0,0)
mixed<-c(1,0)

#Major Gov vs. Minor Opp
data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=mixed,gov=mixed, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]


#Major Gov vs. Minor Gov
data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=mixed,gov=gov, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],conf.level=.95,num.draws=1e5)[1:4]

#Major Opp vs. Minor Opp
data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=mixed,gov=nogov, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],conf.level=.95,num.draws=1e5)[1:4]

#Major Gov vs Major Opp
data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=majorparty,gov=mixed, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]

#Minor Gov vs Minor Opp
data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=minorparty,gov=mixed, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],conf.level=.95,num.draws=1e5)[1:4]

###Seat Share

data.x.mat<-data.frame(Time=Time,seatdiff=seatdiff.10,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=minorparty,gov=nogov, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],conf.level=.95,num.draws=1e5)[1:4]


data.x.mat<-data.frame(Time=Time,seatdiff=seatdifften,groups=groups4,perfemmp2=perfemmp2,seatshare=seatshare,Country=Country,majorparty=minorparty,gov=nogov, FoundedAfter1980= FoundedAfter1980)
x.mat<-build.x.mat(model=model1,data.x.mat)

db.2(model=model1,x.1=x.mat[1,],x.2=x.mat[2,],conf.level=.95,num.draws=1e5)[1:4]


#######

#######
#####
##GENERATING PREDICTED VALUES FOR MODEL 2-PARTY LEADER DEATH
######

LeaderYear<-c(5,5)
nogov<-c(0,0)
opp<-c(0,0)
gov<-c(1,1)
gov.opp<-c(1,0)
minorparty<-c(0,0)
majorparty<-c(1,1)
seatdiff.mean<-c(-0.004717,-0.004717 )
seatdiff.median<-c(0,0)
seatdiff.20<-c(.2,.2) 
seatdiff.15<-c(0.15,.15) 
seatdiff.neg15<-c(-0.15,-.15) 
seatdiff.neg10<-c(-0.10,-.10) 
seatdiff.neg5<-c(-0.5,-.5) 
#seatdiff<-c(.38,-.57) #max and min
#seatdiff<-c(0.0261,-0.0336) #quartile
#seatdiff<-c(0.07395498,-0.07565408) #decile
seatdiff.10<-c(0.10,.10) 
seatdiff.neg10<-c(-0.10,-0.10) 
seatdiff.neg20<-c(-0.20,-0.20) 
seatdiff.neg12<-c(-0.12,-0.12) 
seatdiff.12<-c(0.12,0.12) 
#seatdiff<-c(-0.002417,-0.002417)
#groups<-c(40,60)
groups<-c(30,30)
#groups<-factor(groups,levels=levels(leader$groups))
perfemmp2<-c(20.84,20.84)
#perfemmp2<-c(16.220, 16.220)
#seatshare<-c(.24,.24)
seatshare<-c(0.12570,0.12570)
age<-c(51,51)
Country="Denmark"
#Country=c("Austria","Sweden")
Female<-c(0,1)
Female2<-c(1,0)
Male<-c(0,0)
died<-c(0,0)
partydeath<-c(0,0)
FoundedAfter1980<-c(0,0)
mixed<-c(0,1)


#Neg 15

data.x.mat<-data.frame(LeaderYear=LeaderYear,seatdiff=seatdiff.neg15,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Female,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980= FoundedAfter1980)

x.mat<-build.x.mat(model= deathmodellinear,data.x.mat)

db.2(model=deathmodellinear,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]

#Neg 10

data.x.mat<-data.frame(LeaderYear=LeaderYear,seatdiff=seatdiff.neg10,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Female,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980= FoundedAfter1980)

x.mat<-build.x.mat(model= deathmodellinear,data.x.mat)

db.2(model=deathmodellinear,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]


#Neg 5

data.x.mat<-data.frame(LeaderYear=LeaderYear,seatdiff=seatdiff.neg5,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Female,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980= FoundedAfter1980)

x.mat<-build.x.mat(model= deathmodellinear,data.x.mat)

db.2(model=deathmodellinear,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]


#Pos 15

data.x.mat<-data.frame(LeaderYear=LeaderYear,seatdiff=seatdiff.15,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Female2,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980= FoundedAfter1980)

x.mat<-build.x.mat(model= deathmodellinear,data.x.mat)

db.2(model=deathmodellinear,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]

#Pos 10

data.x.mat<-data.frame(LeaderYear=LeaderYear,seatdiff=seatdiff.10,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Female2,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980= FoundedAfter1980)

x.mat<-build.x.mat(model= deathmodellinear,data.x.mat)

db.2(model=deathmodellinear,x.1=x.mat[1,],x.2=x.mat[2,],num.draws=1e5)[1:4]


##########
##Building Figure 3 
##########
LeaderYear<-5
nogov=0
gov=1
minorparty=0
majorparty=1
seatdiff.mean=-0.002417
seatdiff.median=0.00000
seatdiff.min=-0.566100
seatdiff.max= 0.383000
seatdiff.1=-0.063 
seatdiff.9=0.0659
groups=30
perfemmp2=19.720
perfemmp2.first= 5.359
perfemmp2.third= 31.310 
age=51
Country="Denmark"
femled=1
menled=0
died=0
partydeath=0
seatshare=0.12570
old=0

sd.vec=seq(-.15,.15,by=.01)

Mixed<-c(1)
x.mat.bb<-data.frame(LeaderYear=LeaderYear,seatdiff=sd.vec,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Mixed,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980=old)
x.mat.bb<-build.x.mat(deathmodellinear,x.mat.bb)

x.bb<-apply(x.mat.bb,1,db.2,model=deathmodellinear,x.2=x.mat.bb[sd.vec==0,],unlink=function(x){exp(x)},num.draws=100000)
x.bb<-matrix(unlist(x.bb),nrow=7)
plot(x.bb[3,]~I(sd.vec*100),type="l",col="black",lty=1,lwd=2, ylim=c(min(x.bb[4:5,]),max(x.bb[4:5,])),xlab="Percentage Difference in Seat Share",ylab="Cox Hazard Ratio")
points(x.bb[6,]~I(sd.vec*100),type="l",col="gray55",lty=1)
points(x.bb[7,]~I(sd.vec*100),type="l",col="gray55",lty=1)

Mixed<-c(0)
x.mat.bb<-data.frame(LeaderYear=LeaderYear,seatdiff=sd.vec,groups=groups,perfemmp2=perfemmp2,Country=Country,majorparty=minorparty,gov=nogov,Female=Mixed,leaderage=age,died2=died,Party_Death=partydeath,seatshare=seatshare,FoundedAfter1980=old)
x.mat.bb<-build.x.mat(deathmodellinear,x.mat.bb)

x.bb<-apply(x.mat.bb,1,db.2,model=deathmodellinear,x.2=x.mat.bb[sd.vec==0,],unlink=function(x){exp(x)},num.draws=100000)
x.bb<-matrix(unlist(x.bb),nrow=7)
lines(x.bb[3,]~I(sd.vec*100),col="black",lty=2,lwd=2, )
points(x.bb[6,]~I(sd.vec*100),type="l",col="gray55",lty=2)
points(x.bb[7,]~I(sd.vec*100),type="l",col="gray55",lty=2)


legend(5, 7.5, c("Female Leader", "Male Leader"), 
                col=c("black","black"), lty=c(1,2),lwd=c(2,2))
        
        
######
###THINKING ABOUT TIME
######


##MODEL 1 (TIME UNTIL FEMALE LEADER)

summary(model1.2<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)+I((Time-20)^2)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

summary(model1.3<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)+I((Time-20)^2)+I((Time-20)^3)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

summary(model1.4<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)+I((Time-20)^2)+I((Time-20)^3)+I((Time-20)^4)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

summary(model1.5<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)+I((Time-20)^2)+I((Time-20)^3)+I((Time-20)^4)+I((Time-20)^5)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))


BIC(model1)
BIC(model1.2)
BIC(model1.3)
BIC(model1.4)
BIC(model1.5)

AIC(model1)
AIC(model1.2)
AIC(model1.3)
AIC(model1.4)
AIC(model1.5)

####MODEL 2 (PARTY LEADER DEATH)


summary(deathmodellinear.2<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)+I((LeaderYearFirst-5)^2)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

summary(deathmodellinear.3<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)+I((LeaderYearFirst-5)^2)+I((LeaderYearFirst-5)^3)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

summary(deathmodellinear.4<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)+I((LeaderYearFirst-5)^2)+I((LeaderYearFirst-5)^3)+I((LeaderYearFirst-5)^4)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))


summary(deathmodellinear.5<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)+I((LeaderYearFirst-5)^2)+I((LeaderYearFirst-5)^3)+I((LeaderYearFirst-5)^4)+I((LeaderYearFirst-5)^5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

BIC(deathmodellinear)
BIC(deathmodellinear.2)
BIC(deathmodellinear.3)
BIC(deathmodellinear.4)
BIC(deathmodellinear.5)

AIC(deathmodellinear)
AIC(deathmodellinear.2)
AIC(deathmodellinear.3)
AIC(deathmodellinear.4)
AIC(deathmodellinear.5)

        
                
                
##############
###ROBUSTNESS CHECKS FOR SI
###########


##Discrete Time Duration Model of Women’s Entrance into Party Leadership (With Majoritarian Control Variable)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+majoritarian
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

library(xtable)
xtable(summary(model1),digits=3)

# Discrete Time Duration Model of Party Leaders’ Exit from Office (With Majoritarian Control Variable)


summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+majoritarian
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)

###Discrete Time Duration Model of Women’s Entrance into Party Leadership (With Fem. Mandate Parties)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1980 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0,],family=binomial(link=cloglog),x=TRUE))

summary(model1)
xtable(summary(model1),digits=3)


#Discrete Time Duration Model of Party Leaders’ Exit from Office (Without Fem. Mandate Parties)

summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)

#Discrete Time Duration Model of Women’s Entrance into Party Leadership (Founding Date: 1965)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1965 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(model1),digits=3)

#Discrete Time Duration Model of Party Leaders’ Exit from Office (Founding Date: 1965)

summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1965
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)


#Discrete Time Duration Model of Women’s Entrance into Party Leadership (Founding Date: 1970)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1970 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(model1),digits=3)


# Discrete Time Duration Model of Party Leaders’ Exit from Office (Founding Date: 1970)

summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1970
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)

#Discrete Time Duration Model of Women’s Entrance into Party Leadership (Founding Date: 1975)

summary(model1<-glm(FFL.new~1
+gov*majorparty
+I(seatdiff*100)
+seatshare
+perfemmp2
+ FoundedAfter1975 
+I(Time-20)
+as.factor(groups)
+Country
,data=leader[leader$AFFL.new==0&leader$fem.mandate!=1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(model1),digits=3)

#Discrete Time Duration Model of Party Leaders’ Exit from Office (Founding Date: 1975)

summary(deathmodellinear<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1975
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[leader$LeaderInclude.wC.wPD==1,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear),digits=3)

##Analysis of Leaders' Tenure in Office (First Female Leaders Only)

ind.keep.DML<-(leader$Female==0 | leader$FemLeaderNum==1) & leader$LeaderInclude.wC.wPD==1

summary(deathmodellinear.ffo<-glm(LeaderDeath~
1
+I(seatdiff*100)*Female
+gov*majorparty
+seatshare
+FoundedAfter1980
+perfemmp2
+leaderage
+died2
+Party_Death
+I(LeaderYear-5)
+as.factor(groups)
+ Country
,data=leader[ind.keep.DML,],family=binomial(link=cloglog),x=TRUE))

xtable(summary(deathmodellinear.ffo),digits=3)


####
###DESCRIPTIVE STATS
#####

###NOTE: The three parties with dual male-female leaders are included twice in the dataset (once for male leader, once for female)
##These three parties should not be doubly counted in descriptive stats


#Total Leaders
leader$countrypartyid[which(leader$NewLeader==1&leader$LeaderInclude.wC.wPD == 1)]

#Female Leaders
summary((xtabs(~Party+ Female ,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

leader$countrypartyid[which(leader$NewLeader==1&leader$Female==1&leader$LeaderInclude.wC.wPD == 1)]

leader$countrypartyid[which(leader$FFL.new==1)]

(xtabs(~Female,data=leader[rownames(model.matrix(deathmodellinear)),]))

#Change in Seat Share

round(summary(leader[rownames(model.matrix(deathmodellinear)),"seatdiff"])*100,3)
round(sd(leader[rownames(model.matrix(deathmodellinear)),"seatdiff"])*100,3)
round(IQR(leader[rownames(model.matrix(deathmodellinear)),"seatdiff"])*100,3)


#Size

round(summary(leader[rownames(model.matrix(deathmodellinear)),"seatshare"])*100,3)
round(sd(leader[rownames(model.matrix(deathmodellinear)),"seatshare"])*100,3)
round(IQR(leader[rownames(model.matrix(deathmodellinear)),"seatshare"])*100,3)


#Per Fem Mp

round(summary(leader[rownames(model.matrix(deathmodellinear)),"perfemmp2"]),3)
round(sd(leader[rownames(model.matrix(deathmodellinear)),"perfemmp2"]),3)
round(IQR(leader[rownames(model.matrix(deathmodellinear)),"perfemmp2"]),3)

#Leader Age

(summary(leader[rownames(model.matrix(deathmodellinear)),"leaderage"]))
(sd(leader[rownames(model.matrix(deathmodellinear)),"leaderage"]))
(IQR(leader[rownames(model.matrix(deathmodellinear)),"leaderage"]))



#Country
apply(xtabs(~Party+Country,data=leader[rownames(model.matrix(deathmodellinear)),])>0,2,sum)

(xtabs(~Country,data=leader[rownames(model.matrix(deathmodellinear)),]))


#Party Family
apply(xtabs(~Party+groups,data=leader[rownames(model.matrix(deathmodellinear)),])>0,2,sum)

(xtabs(~groups,data=leader[rownames(model.matrix(deathmodellinear)),]))



#Government
summary((xtabs(~Party+gov,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

(xtabs(~gov,data=leader[rownames(model.matrix(deathmodellinear)),]))

#Major Party
summary((xtabs(~Party+majorparty ,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

(xtabs(~majorparty ,data=leader[rownames(model.matrix(deathmodellinear)),]))

#FoundedAfter1980
summary((xtabs(~Party+ FoundedAfter1980 ,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

(xtabs(~ FoundedAfter1980 ,data=leader[rownames(model.matrix(deathmodellinear)),]))

#Party_Death
summary((xtabs(~Party+ Party_Death ,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

(xtabs(~ Party_Death ,data=leader[rownames(model.matrix(deathmodellinear)),]))

#Died
summary((xtabs(~Party+ died2 ,data=leader[rownames(model.matrix(deathmodellinear)),])>0)[,2])

(xtabs(~ died2 ,data=leader[rownames(model.matrix(deathmodellinear)),]))



